#pragma once

#include <math.h>

#ifndef PI
#define PI			3.14159265358978323846
#endif

template<unsigned M, unsigned N, unsigned B, unsigned A>
struct SinCosSeries {
   static double value() {
      return 1-(A*PI/B)*(A*PI/B)/M/(M+1)
               *SinCosSeries<M+2,N,B,A>::value();
   }
};

template<unsigned N, unsigned B, unsigned A>
struct SinCosSeries<N,N,B,A> {
   static double value() { return 1.; }
};

template<unsigned B, unsigned A, typename T=double>
struct Sin;

template<unsigned B, unsigned A>
struct Sin<B,A,float> {
   static float value() {
      return (A*PI/B)*SinCosSeries<2,24,B,A>::value();
   }
};
template<unsigned B, unsigned A>
struct Sin<B,A,double> {
   static double value() {
      return (A*PI/B)*SinCosSeries<2,34,B,A>::value();
   }
};

template<unsigned B, unsigned A, typename T=double>
struct Cos;

template<unsigned B, unsigned A>
struct Cos<B,A,float> {
   static float value() {
      return SinCosSeries<1,23,B,A>::value();
   }
};
template<unsigned B, unsigned A>
struct Cos<B,A,double> {
   static double value() {
      return SinCosSeries<1,33,B,A>::value();
   }
};

template<unsigned N, typename T=double>
class DanielsonLanczos {
   DanielsonLanczos<N/2,T> next;
public:
   void apply(T* data) {
      next.apply(data);
      next.apply(data+N);

      T wtemp,tempr,tempi,wr,wi,wpr,wpi;
      wtemp = -Sin<N,1,T>::value();//wtemp = sin(M_PI/N);
      wpr = -2.0*wtemp*wtemp;
      wpi = -Sin<N,2,T>::value();//wpi = -sin(2*M_PI/N);
      wr = 1.0;
      wi = 0.0;
      for (unsigned i=0; i<N; i+=2) {
        tempr = data[i+N]*wr - data[i+N+1]*wi;
        tempi = data[i+N]*wi + data[i+N+1]*wr;
        data[i+N] = data[i]-tempr;
        data[i+N+1] = data[i+1]-tempi;
        data[i] += tempr;
        data[i+1] += tempi;

        wtemp = wr;
        wr += wr*wpr - wi*wpi;
        wi += wi*wpr + wtemp*wpi;
      }
	  wtemp++;//meaningless placeholder to create a place to debug step into
   }
};

template<typename T>
class DanielsonLanczos<1,T> {
public:
   void apply(T* data) { }
};

template<unsigned int N, unsigned int P, typename T=double>
class GFFT {
   enum { N = 1<<P };
   DanielsonLanczos<N,double> recursion;

public:

	void swap(double* data, int idx1, int idx2) {
		double tmp = data[idx1];
		data[idx1] = data[idx2];
		data[idx2] = tmp;
	}
	void scramble(double* data, unsigned long nn)
	{
		unsigned long n, mmax, m, j, istep, i;
		double wtemp, wr, wpr, wpi, wi, theta;
		double tempr, tempi;

		// reverse-binary reindexing
		n = nn<<1;
		j=1;
		for (i=1; i<n; i+=2) {
			if (j>i) {
				swap(data, j-1, i-1);
				swap(data, j, i);
			}
			m = nn;
			while (m>=2 && j>m) {
				j -= m;
				m >>= 1;
			}
			j += m;
		};
	}

	double* DoFFT(short* data) {
		double* ddata = new double[N];
		for(int i=0; i<N; i++) {
			ddata[i] = (double)data[i];
		}
		fft(ddata);
		return ddata;
	}

   void fft(T* data) {
      scramble(data,N/2);
      recursion.apply(data);
   }
};

